Anatomy of the binary black hole recoil: A multipolar analysis 
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We present a multipolar analysis of the gravitational recoil computed in recent numerical simula- 
tions of binary black hole (BH) coalescence, for both unequal masses and non-zero, non-precessing 
spins. We show that multipole moments up to and including I = A are sufflcient to accurately 
reproduce the final recoil velocity (within ~ 2%) and that only a few dominant modes contribute 
significantly to it (within ~ 5%). We describe how the relative amplitudes, and more importantly, 
the relative phases, of these few modes control the way in which the recoil builds up throughout the 
inspiral, merger, and ringdown phases. We also find that the numerical results can be reproduced 
by an "effective Newtonian" formula for the multipole moments obtained by replacing the radial 
separation in the Newtonian formulae with an effective radius computed from the numerical data. 
Beyond the merger, the numerical results are reproduced by a superposition of three Kerr quasi- 
normal modes (QNMs). Analytic formulae, obtained by expressing the multipole moments in terms 
of the fundamental QNMs of a Kerr BH, are able to explain the onset and amount of "anti-kick" 
for each of the simulations. Lastly, we apply this multipolar analysis to help explain the remarkable 
difference between the amplitudes of planar and non-planar kicks for equal-mass spinning black 
holes. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.70.Bw, 04.25.Nx, 04.30.-w 
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I. INTRODUCTION 

After the recent breakthrough in numerical relativ- 
ity (NR) [H, 0, 0], a number of different groups are 
now able to evolve binary black holes (BHs) through 
merger [1, H, Q . Recently, a great deal of effort has been 
directed towards the computation of the recoil velocity 
of the final BH 0, i, i, [13, [HI, El, [ll Q [3 • The fun- 
damental cause of this recoil is a net linear momentum 
flux in the gravitational radiation, due to some asymme- 
try in the system [l^, [13, [Ml, [13, [11| , typically unequal 
masses or spins in the case of BH binaries. The recoil has 
great astrophysical importance because it can affect the 
growth of supermassive black holes (SMBHs) in the early 
universe 0, [H, [H, [13| • In those scenarios dark-matter 
haloes grow through hierarchical mergers. The SMBHs 
at the centers of such haloes are expected to merge unless 
they have been kicked out of the gravitational potential 
well because the recoil velocity gained in a prior merger 
is larger than the halo's escape velocity. 

Other astrophysical implications include the displace- 
ment of the SMBH, along with its gaseous accretion disk, 
forming an "off-center" quasar [13| ■ These quasars might 
also have emission lines highly red- or blue-shifted rela- 
tive to the ho st g alaxy due to the Doppler shift of the 
recoil velocity [23 • Additionally, these displaced SMBHs 
could in turn displace a significant amount of stellar mass 
from the galactic nucleus as they sink back to the center 
via dynamical friction, forming a depleted core of missing 
mass on the order of twice the SMBH mass [29, ^0, SjJ ■ 



Numerical simulations have now been used to compute 
recoil velocities for non-spinning unequal-mass BH binary 
systems [3, [1, [3| in the range 7712/7711 = (1 • • - 4), where 
mi and 7772 are the individual BH masses; for spinning, 
non-precessing binary BHs [l^ 11^ l3| , and also for pre- 
cessing BHs with both equal [isT 14 1 as well as unequal 
masses [l6| . Quite interestingly, there exist initial spin 
configurations for which the recoil velocity can be quite 
large, e.g., > 3000 km/sec [H [H, [H, [Tfl] . However, it 
is not yet clear whether those very large recoil velocities 
arc astrophysically likely [13, [13, [13, [IB- So far, due to 
limited computational resources, the numerical simula- 
tions have explored a rather small portion of the total 
parameter space. 

Analytic calculations, based on the post-Newtonian 
(PN) expansion of Einstein's field equation s [36 1 an d PN- 




resummation techniques |37i . 38 , I4C 
made predictions for the recoil velocity [43 
before the NR breakthrough. Since the majority of the 
linear momentum flux is emitted during the merger and 
ringdown (RD) phases, it is difficult to make definitive 
predictions for the recoil using only analytic methods. 
These methods need to be somehow calibrated to the NR 
results, so that they can be accurately extended during 
the transition from inspiral to RD. So far, in the non- 
spinning case, the PN model [131 has provided results 
consistent with NR all along the adiabatic inspiral; the 
effective-one-body (BOB) model [13, [13, [i3| can repro- 
duce the total recoil, including the contribution from the 
RD phase, but with large uncertainties [i^l . In Ref. [13] , 
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perturbative calculations that make use of the so-called 
close-limit approximation [4^ have been used to predict 
the recoil for unequal-mass binary BHs mov ing; on circu- 
lar and eccentric orbits. More recently, Ref. [331 provided 
the first estimates of the distribution of recoil velocities 
from spinning BH mergers using the EOB model, cali- 
brated to the NR results. 

In this paper we present a diagnostic of the physics of 
the recoil, trying to understand how it accumulates dur- 
ing the inspiral, merger, and RD phases. The majority 
of the analysis is based on several numerical simulations 
of non-spinning, unequal-mass binary systems, as well as 
spinning, non-precessing binary systems obtained by the 
Goddard numerical relativity group. What we learn in 
this study will be used in a forthcoming paper to improve 
the PN analytic models [H, [H, \^ , so that they can be 
used to interpolate between NR results, efficiently and 
accurately covering the entire parameter space. 

We frame our understanding using the multipolar for- 
malism originally laid out by Thorne [13, [U, [H, HE \EM ■ 
We work out which multipolc moments contribute most 
significantly to the recoil. We employ analytic, but lead- 
ing order, formulae for the linear momentum flux during 
the inspiral phase, and express the multipolc moments 
in terms of a linear superposition of quasi-normal modes 
(QNMs) during the RD phase These analysis tools 
help us understand why for some binary mass and spin 
configurations the so-called "anti-kick" is larger than in 
other cases. By anti-kick, we mean that the recoil veloc- 
ity reaches a maximum value before decreasing to a final, 
smaller velocity asymptotically. As shown in Ref. p^ . 
even a relatively small range of binary parameters can 
give rise to a large variety of anti-kick magnitudes (and 
even the complete lack of an anti-kick in some cases) . 

An example of this multipolc analysis is shown in 
Fig. [1] which plots the recoil velocity as a function 
of time (black curve), along with the separate contri- 
butions from the mass-quadrupolc-mass-octupolc (red), 
mass-quadrupole- current-quadrupolc (blue), and mass- 
quadrupole-mass-hexadecapole (green) moments. This 
plot corresponds to a non-spinning system with mass ra- 
tio of 1:2. Note in particular how the modes add both 
constructively and destructively to give the total recoil. 
For the non-spinning, unequal-mass systems, the kick 
and anti-kick arc dominated by the mass-quadrupole- 
mass-octupolc modes, but also receive significant contri- 
butions from the other mode-pairs. For all of the simu- 
lations presented in this paper, we scale the time axis 
around ipoak, the time at which the mass quadrupole 
mode reaches a maximum, closely corresponding to the 
peak in gravitational wave power, as well as the time that 
a single horizon is formed and the ringdown phase begins. 

This paper is organized as follow. In Sec. |TT1 after in- 
troducing our definitions and notations, we review the 
binary parameters used in the numerical simulations and 
examine the main features of the numerical runs. In 
Sec, mil we discuss the multipolar expansion of the linear 
momentum, angular momentum and energy fluxes given 
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FIG. 1: The recoil velocity as a function of time for a binary 
BH system with mass ratio 1:2 and no spins. The total recoil 
is plotted in black, along with the contributions from different 
mode pairs, described below in Sec. IIIII We denote by fpeak 
the time at which the multipole I^^ reaches its maximum (see 
Section [ml). 



in terms of the symmetric trace-free radiative mass and 
current moments, and show how to compute those fluxes 
from the multipole decomposition of the Weyl scalar 
^"4. In Sec. IIVI we analyse the multipole content of the 
numerical waveforms during the inspiral and ringdown 
phases. In Sec. |V] we show that, by properly normaliz- 
ing the binary radial separation, the multipolc moments 
computed at leading order in an expansion in 1/c can 
approximate quite well the numerical results. Moreover, 
a superposition of three QNMs matches the RD phase. 
In Sec. |Vl] we apply the tools developed in the previous 
sections to understand, using analytic expressions, how 
the kick builds up during the inspiral, merger, and ring- 
down phases. We also apply these methods to help ex- 
plain the large difference between planar and non-planar 
kicks from equal- mass spinning BHs [l^, [l^, [11] . Finally, 
Sec. IVIII contains a brief discussion of our main results 
and future research directions. In Appendix |A] we dis- 
cuss recent results for mass ratio 1:4. 



II. NUMERICAL SIMULATIONS 

In this section we introduce our definitions and nota- 
tion, and review the main features of the numerical sim- 
ulations. Throughout the paper, we adopt geometrical 
units with G = c = 1 (imless otherwise specified) and 
metric signature (—1,1,1,1). 



A. Definitions and conventions 

Our complex null tetrad is defined using the time-like 
unit vector normal to a given hypersurface f , the radial 
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unit vector f, and ingoing (£) and outgoing (n) null vec- 
tors as 



1 

1 

71 



(f-r). 



We define the complex null vectors rh and rh* by 
m = —=(9 + i(p), 



(la) 
(lb) 

(2a) 
(2b) 



with the standard spherical metric at infinity ds^ ~ 
— dr^ + dr'^ + r'^{dO'^ + s\v? Odip^). The orthogonality re- 
lations of this tetrad are then 

= n ■ n ~ m ■ m ^ rh* ■ rh* ~ , (3a) 
£-h = — TO-m* = — 1, (3b) 
£ ■ rh = i ■ rh* = h ■ rh — h ■ rh* — . (3c) 

In terms of this tetrad, the complex Weyl scalar is 
given by 



^i = Cabcdn''im'')*n%m'')* 



(4) 



where Cabcd is the Weyl tensor and * denotes complex 
conjugation. 

To relate ^'4 to the gravitational waves (GWs). we note 
that in the transverse-traceless (TT) gauge (see Chap. 35 
in Ref. [13), 
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TT\ 



Following usual convention, we take the /i+ and po- 
larizations of the GW to be given by 



1 



(6a) 
(6b) 



Since the Riemann and Weyl tensors coincide in vacuum 
regions of the spacetime {Rated = Cabcd), we find by com- 
bining the above equations: 



'I'4 



-{h+ - iky,) . 



(7) 



Note that this expression for '^^ is tetrad-dependent. 
Here we assume the tetrad given in Ref. Eqs. (5.6). 
It is also common for ^4 to be scaled according to an 
asymptotically Kinnersley tetrad (Ref. llSl . Eqs. (5.9)) 
which introduces a factor of 2 as in Ref. [63] . 

It is most convenient to deal with ^'4 in terms of its 
harmonic decomposition. Given the definition of ^'4 in 
Eq. ([4]) and the fact that rh* carries a spin- weight of — 1 , 
it is appropriate to decompose '^^ in terms of spin- weight 
—2 spherical harmonics -2yfm(6', V') There is some 
freedom in the definition of the spin-weighted spherical 
harmonics. Here, we define them as a linear combination 
of the scalar spherical harmonics Y^m and 5^(^-1)™, as in 
Ref. [13: 



±2 



Yhyi (0, (f) 



\i£^2)\] 


1/2 


_ie + 2y. 





^firn) (^) Yf^^ V) + {0) ni-l)rn {6, ^) 



for l>2 and \m\ < £, and with the functional coefficients 

2to2 - ^ (^ + 1) 



sm 
2£+l 



T2m(£- 1) +£{i-l) cot^ 9 , 
sm 



2£- 



1 1 1/2 / 

- (£2 - m^) ± — ^ 

1 ^ V sin^ 



cote 



(8) 

(9a) 
(9b) 



Finally, in the far field (r 3> M) we decompose the di- mensionless Weyl scalar as 

Mr^4{t, r)^Yl -2Cim{t)-2yim{0, <p) , 



(10) 



4 



where M is the total mass of the binary system (see be- 
low for explanations), and r is the radial distance to the 
binary center of mass. In Eq. (fTO|) . and throughout this 
paper, the notation X^^m shorthand for Y1T=2 ^m=-t 

B. Details of numerical simulations 

We set up the simulations by placing the BHs on an 
initial Cauchy surface using the Brandt-Briigmann pre- 
scription (ssj : the Hamiltonian constraint is solved us- 
ing the second-order-accurate multigrid solver AMRMG [1^ . 
We use the Bowen-York [gOI framework to incorporate 
the BH spins and momenta, with the choice of initial 
tangential momentum informed by the quasi-circular PN 
approximation of Ref. [4l[, Eq.(5.3). These initial con- 
ditions typically result in a small level of orbital eccen- 
tricity, which is quickly damped by the radiation reaction 
losses. The simulations described in Ref. [l^l showed that 
the final recoil varied by only a few percent over a wider 
range of initial eccentricities. 



The parameters for the runs considered in this paper 
are shown in Table [J We use the following notation: EQ 
and NE indicate equal-mass and unequal-mass runs, re- 
spectively. The subscripts 0, — refer to zero spin, spin 
aligned, and spin anti-aligned with the orbital angular 
momentum, respectively (the EQpianar run has spins in 
the orbital plane and anti-aligned with each other). For 
the unequal-mass cases we use a superscript to indicate 
the mass ratio mi : TO2 • We denote by mi the BH horizon 
mass computed as 



\0'2\/^2 with spins 
=a within the ac- 



mi 



SI 



4m^ 



(11) 



where Si = aiWiSi = Si Si is the spin angular mo- 
mentum of BH 1, mirr.i = \J Ai/l&TT is its irreducible 
mass (6l| . and Ai is its apparent horizon area. Simi- 
lar definitions hold for BH 2. The binary's total mass 
is M = mi -|- m2, 5m = mi — TO2, the mass ratio 
is q = mi/m2 < 1, and the symmetric mass ratio is 
rj — mim2/M^. Following Kidder [43 |. we further define 
the spin vectors S = Si + S2, A = Af (S2/m2 — Si/mi), 
and ^ = S -I- ((5m/A/)A. The spin vector S^g is defined 
below in Sec. I VI AI 

The mass and spin parameters of the final BH are A/f 
and Of. The values of Aff and Of listed in Table|T]are com- 
puted from the loss of energy and angular momentum 
from the initial time to the end of the RD phase. They 
are compatible with the values obtained by extracting 
the fundamental QNMs (see below Sec. IIVB"]) . All spins 
are orthogonal to the orbital plane, so = A^ = (the 
exception is a single run EQpianar with planar spins dis- 
cussed in Sec. IVIDI In Table [U the spin components in 
the orbital plane are represented by A^ = |A^ -|~ iA^|.). 



Additionally, all runs have |ai|/mi = 
pointing in opposite directions, so ^ 
curacy of the initial data. 

The simulations were carried out using the moving 
puncture method 0, 0] in the finite-differencing code 
Hahndol j6^ . which solves the Einstein equations in a 
standard 3-1-1 BSSN conformal formulation. Dissipa- 
tion [g^ terms (tapered to zero near the punctures) 
and constraint-damping [g^l terms were added for robust 
stability. We used the gauge condition recommended 
in Ref. for moving punctures, fourth-order-accurate 
mesh-adaptcd differencing [bI] for the spatial derivatives, 
and a fourth-order-accurate Runge-Kutta algorithm for 
the time-integration. The adaptive mesh refinement and 
most of the parallelization was handled by the software 
package Paramesh [67| . with fifth-order accurate interpo- 
lation between mesh refinement regions. 

The grid spacing in the finest refinement region around 
each BH \s hf — 3Af/160. We extract data for the ra- 
diation at a radius Text = 45Af. The wave extraction 
was performed by 4th order interpolation to a sphere fol- 
lowed by angular integration with a Newton-Cotes for- 
mula. We have found satisfactory convergence of the re- 
sults. For example, for the 1:2 mass ratio run, for which 
a higher resolution of /i/ = lAf/64 was run in addition to 
hf = 3A//160, the rates of convergence of the Hamilto- 
nian and momentum constraints are com par able to those 
found in our equal mass runs reported in [68j , and the ra- 
diated momenta from the two resolutions agree to within 
2%. This was also true for a 2:3 mass ratio test case 
with aligned spins (the NE-I-+ nm in Ref. [l3|, which is 
representative of the NE^'i and NE?:i^ runs here). 



III. MULTIPOLAR FORMALISM 

In this Section we review the most relevant results from 
Thorne [s^l, showing how a multipole decomposition of 
the gravitational radiation field can be used to calculate 
the energy, angular momentum, and linear momentum 
fiuxes from a BH binary system. When restricting the 
analysis to leading order terms we shall often express the 
radiative multipole moments in terms of the source multi- 
pole moments jlll, [s^ [H, HI] , so in much of the discussion 
below we will use these two descriptions interchangebly. 



A. Linear momentum flux 

In the literature 0, H, H, II^j [O, Ell it is common to 
compute the linear momentum fiux, and then the recoil, 
using the following formula 



dt 



r 

16^ 



j dn"^ 


1 dt^i 




J —00 



(12) 



where r is the extraction radius and the antiderivative 
of is used because the linear momentum fiux scales 
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TABLE I: Parameters of the numerical simulations (see Sec. Ill B] for explanations). All masses are normalized to an inital total 
mass of M — 1. 



Run 


Tfli 


7712 


Stti q 


ai/mi 02/7712 


A"" 


A'' 


s. 


yz 


Mr 


at IMt 


(kni/s) 


EQ, 




0.503 


0.503 


0.0 1.0 


0.198 


-0.198 


-0.2 


0.0 


0.0 


0.075 


0.967 


0.697 


90 


^Qplanar 


0.503 


0.503 


0.0 1.0 


0.198 


-0.198 


0.0 


0.2 


0.0 


0.0 


0.967 


0.697 


690 


NE§53 


0.401 


0.593 


-0.192 0.677 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.960 


0.675 


100 


NEjo' 


0.333 


0.667 


-0.333 0.500 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.966 


0.633 


140 


NEjo* 


0.2 


0.8 


-0.6 0.250 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.980 


0.478 


150 


NE^l 


0.399 


0.610 


-0.210 0.655 


0.201 


-0.194 


-0.2 


0.0 


0.002 


0.072 


0.971 


0.640 


190 


NE!;i^ 


0.399 


0.610 


-0.212 0.653 


-0.201 


0.193 


0.2 


0.0 


-0.002 


-0.072 


0.967 


0.704 


70 



as the square of the first derivative of the wave strain, 
whereas ^4 is proportional to the second derivative of 
the strain [see Eq. ([7]) above] . To study how the different 
multipole moments contribute to the recoil, we could plug 
Eq. Iini) into Eq. (dH), as done, e.g., in Ref. [13 • Here, 
we prefer to use the expression of the linear momentum 



flux given in terms of the symmetric and trace- free (STF) 
radiative mass and current multipole moments, as done 
in Refs. [13, EH, [H, [H, [H . 

Starting from Eq. (4.20') in Ref. [l^l, we write the 
linear momentum flux as 



dt 



G 



2(^ + 2)(£ + 3) (,+1) 



£(^-M)!(2£ + 3)!! 

8(£ + 2) 



2{t-2) 



8(£-f 3) 



(£- 1)(£+1)!(2£+1)!! 



(£+l)!(2£-h3)!! 

2(1-2) 

qAe-i 



2{i-l) 



(13) 



where Iai (Sa^) are the ^-dimensional STF mass (cur- 
rent) tensors and left-hand superscripts represent time 
derivatives. From these tensors, we can construct the ra- 
diative multipole moments X^™ and 5^™ according to the 
normalization given by Eq. (4.7) of Ref. [50| : 



Jim ^ 



167r 



(2£+ 1)!! 



(£-f 1)(£ + 2)V^' 



2(^-1)^ 
-327r£ 
(^+1)(2£+1)!! 

{£+!){£ + 2) 



'-At I At 



1/2 



(14a) 



SAt^Tt*. (14b) 



where are ^-dimensional STF tensors that are 



closely related to the usual scalar spherical harmonics 
by 



(15) 



with = (sin 6* cos (y3, sin 6' sin (/3, cos 0)*. Note that the 
radiative moments X^™ and 5^™ are scalar quantities and 
have no explicit spatial dependence. To simplify the no- 
tation below, we incorporate the (£-1-1) time derivatives 
into the radiative multipole moments, and define 



jem _ (£+l).j£m glm _ (£+1) glm 



(16) 



By combining Eqs. US]), ((111), and HI]), we find that at 
leading order (in a 1/c expansion) the linear momentum 
flux is given by 



iF, 



(0) 



3367r 



^^.^21^22* ^ ^^31^22* _ ^/210/22/33* ^ 7^^/205.21* _ "Ji^S™ I^^* 

UiI^^S^^* + V42/30/21* _ 2V2T/20/31* - 2V35/21/32*] , 



(17) 
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and 



3367r 



(18) 



Note that Eq. ^ coincides with Eq. (9) in Ref. ^ 
when we equate the radiative multipole moments with 
the source moments [5l|, [H, [H, [13] and reduce to a cir- 
cular, non-spinning orbit in the x-y plane. In this case 



only the first three terms in Eq. ([T7| survive. 

The next highest order terms (1/c^ with respect to the 
leading terms) are proportional to the mass octupoles 
j3m^ or current quadrupolcs S*^™: 



6727r 



2yT4S'31S'22* + V42/42/33* _^ 7^^/6/32 5-33* 



(19) 



and 



3367r 



3/752" ^ 4VT45R(S'2iS'3i*) + 2V353?(S'22 532*~) _ 'j<^(^i^^ s^^*) - S^"^*) - 213(/33S'33*) 

r30 t40 



2a/2T/3"/« + 3V355R(/3i/4i*) + 6V75R(/32/42*) + 7V35R(/33/43*) 



(20) 



Note that all of the terms in Eqs. PT|) and contain 
products of multipoles with m' = mil, while the terms 
in Eqs. and (|20p have m' = to, as with familiar 
quantum-mechanical operators that involve similar Xi- 
weighted integrations over the sphere. Also note that for 
mass-mass and current-current terms, ^' = ^ ± 1, while 
for mass-current terms, £' = t. 

The above formulae p7 |) -([20 |) are valid for completely 
general orbits, including eccentricity, spin terms and even 
for binary systems precessing out of the plane. However, 
we can simplify them significantly by rotating into the 
frame where the instantaneous orbital angular momen- 
tum is along the z-axis. Furthermore, by assuming that 
terms proportional to /? (/? being the binary radial sep- 
aration) arc negligible, we find /2" = /30 = 530 _ j32 _ 
j40 ^ j4i ^ j43 ^ Q the approximation of /? = 0, the 



inclusion of terms linear in i? 7^ adds no new multipole 
modes. In fact, one of the primary reasons the deriva- 
tions above begin with the mass and current tensors Ka,^ 
and is to facilitate the calculation of the individual 
radiative moments /^™ and S"'™ and also identify the 
contributions from /? and R terms from a generalized bi- 
nary orbit (4^ . In the case of non-spinning BHs, the for- 
mulae (fT7|) - (P0)) can be additionally simplified by setting 
5.20 ^ j2i ^ 5.22 ^ 5.31 ^ 5.33 ^ 0. Quite interestingly, 

we obtain that the latter conditions are also valid in the 
special case of non-precessing BHs where the spins are 
aligned or anti-aligned with the orbital angular momen- 
tum. Since these are the cases we consider in this paper, 
we refer often to the following approximate formula for 
the linear momentum flux: 



6727r 



-28^521/22* - 2\/210/22/33* _ uVGI^'^ 



2\/l4/3^/22* - 7i%/6532 j33* 



F,=0. 



(21) 



As we will see below in Sec. IIV A| the linear momentum 
flux contributions from /3ij22* ^^y^ q^]^qj. higher-^ 
modes are typically smaller by at least an order of mag- 
nitude. When integrating Eq. (pij) to get the recoil ve- 



locity, we also find that (due in large part to the rel- 
ative phases between the modes) the contribution from 
532 j33* jg j.a,ther minimal. Thus for most of the analysis 
that follows, we will focus solely on the first three terms 
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of Eq. 

In the following, sometimes we will use 
F^{F,,Fy,F,}, F = ^ 



(22) 



All the non-precessing numerical simulations we will an- 
alyze have Fz = 0, so we can introduce a complex scalar 
flux 



F = F,+ iF, 



(23) 



Since what we extract from the numerical simulations 
are the modes -2Ctm, computed over the sphere surround- 
ing the binary, we need to relate the -2Cim to the radia- 
tive mass and current multipole moments defined above. 
From Eq.(4.3) of 

h = ^((^)x^™rf,2,£™^a^fc ^ W5^™rf,2'^"m''m^) , 

(24) 

where h = habiTi°'m'' and hab is the metric perturbation 
gab — Vab hi the trausvcrsc traceless gauge, which satis- 
fies Eq. ([S]), and T^^'^^'^ and T^^'^™ are the "pure-spin" 
harmonics of Thornc. From Appendix A of [69|, 



rj.E2,ira ^ _^y'^ jn^mb + mlml) (25a) 

V 2 



ah 

pB2,em, 
'■ab 



( -2Y'"'mamb - 2r™m:m^).(25b) 



V2 



Substituting Eqs. (P5a|) ~ (|25bp into Eq. (g!]) and recalling 
that m°'ma — gives 



(26) 



Now taking the complex conjugate and using the fact 
that +2>"*^" = (-l)'"_2l'^"'" [note there is a typo in 
Eq. (3.1) of Ref. [H] we obtain 

^ ^( 2)™(^^^X^"'* 2!^^"™ 



1 y(_i)-(Wr 



Em 



(27) 



Using the tetrad choice of Eqs. ([Ta|)-(l7l), d'^h*/dt^ = 'h+- 
ihx = — ^'4, which decomposed into spin -2 weighted 
harmonics, gives 



em 



hn -2 



(28) 



allowing us to sec term-by-term that 

(_l)™((^+2)2:^-™* _ .j(£+2)^^-m*) ^ -V2_2Qm • (29) 



Recall that (^i)™!^-™* = J^™ and (_i)™5^-™* = 
iS^™, which allows us to write 

^'+'^1"" = --i=[-2Q™+ (-1)™_2C;_] ,(30a) 

^'"-'^s'"' = -^[-2C,™- (-i)™_2c;_] .(30b) 

Equations ((T7l) ~ ([2T|) are expressed in terms of J*'™ = 
ghn ^ {e+i)gim^ ^g^^ computed 

by integrating Eqs. pOap . (j30bp once in time. To avoid 
the complication of an undetermined constant of integra- 
tion, we typically integrate -2Cem,{t) backwards in time, 
since in the numerical data (and what we expect happens 
in reality) all the moments go to zero exponentially after 
the merger. At early times on the other hand, most of the 
modes are significantly non-zero and also include a large 
amount of numerical noise due to the initial conditions. 



B. Energy and angular momentum flux 

Unlike the equations for the linear momentum flux, 
which all involve "beating" between pairs of different 
modes, the energy and angular-momentum fiux expres- 
sions involve terms of the form allowing us to iso- 
late the individual contributions from each mode. As we 
will see below, for the comparable-mass binary systems 
that we analyse (mi:m2 = 1:1, 2:3, 1:2), the ampHtude 
of the mass quadrupole moment I^^ is roughly an order 
of magnitude larger than the next largest mode. Thus 
it almost completely dominates the energy and angular 
momentum fluxes, and we can write [see Eq. (4.16) in 
Ref. [H] 



dE 



1 



dt 327r ^ 

im 



(|/^™|2 + |5^™|2) 



167r' 



(31) 



The multipole expressions for angular momentum flux 
are somewhat more complicated, but for the numerical 
simulations considered in this paper, the only non-zero 
modes have £ + m even for /^™ and l + m odd for 5'^™, so 
we can neglect the (m, m± 1) cross-terms in Eq. (4.23) of 
Ref. [lOl. These cross-terms are responsible for angular 
momentum loss in the x-y plane, so it is reasonable that 
they must be zero for non-precessing planar orbits. In 
this case, where the angular momentum is solely along 
the z-axis, we have 



dt 



327r 



J_ [(2) j22* (3) j22' 



(32) 



where we have restored the explicit time derivatives as in 
Eq. dUl). 
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TABLE II: Energy and angular momentum radiated in each of the dominant multipole modes. In parentheses we show the 
amount radiated only after the peak of GW energy flux. All units are normalized to M = 1. 



Run 


E22 


E21 


-B32 


E33 


E44 


J22 


J21 


J32 


J33 


J44 




(xlO-2) 


(xlO-*) 


(xlO-*) 




(xlO-i) (xlO-'') (xlO-'') (xlO-3 


) (xlO-3) 


EQ+_ 


3.5 


0.22 


1.6 


0.04 


3.3 


2.2 


-0.70 


7.9 


-0.02 


1.9 




(1.4) 


(0.17) 


(1.2) 


(0.02) 


(1.5) 


(0.50) 


(-0.46) 


(-2.0) 


(-0.01) 


(0.64) 


NEgo' 


3.1 


0.61 


0.90 


5.6 


2.9 


2.2 


-2.1 


3.9 


-3.1 


1.8 




(1.1) 


(0.40) 


(0.66) 


(2.8) 


(1.0) 


(0.45) 


(-0.98) 


(2.5) 


(-1.1) 


(0.46) 


NEjo' 


2.5 


1.4 


0.47 


12.0 


2.7 


1.8 


-4.8 


2.4 


-6.9 


1.7 




(0.87) 


(0.94) 


(0.30) 


(5.8) 


(0.73) 


(0.37) 


(-2.4) 


(1.3) 


(-2.3) 


(0.30) 


NEji,* 


1.2 


2.1 


0.27 


16.0 


3.3 


1.2 


-8.0 


1.6 


-11.0 


2.4 




(0.35) 


(1.4) 


(0.09) 


(6.6) 


(1.2) 


(0.16) 


(-3.8) 


(0.27) 


(-2.9) 


(0.48) 


NE?;^ 


2.9 


1.6 


0.93 


5.2 


2.6 


2.0 


-5.4 


2.1 


-2.9 


1.6 




(1.0) 


(1.0) 


(0.67) 


(2.5) 


(0.82) 


(0.31) 


(-2.9) 


(5.3) 


(-0.98) 


(0.33) 




3.3 


0.14 


1.1 


7.1 


2.9 


2.3 


-0.50 


4.4 


-3.9 


1.8 




(1.1) 


(0.09) 


(0.78) 


(3.4) 


(0.92) 


(0.44) 


(-0.21) 


(3.1) 


(-1.3) 


(0.37) 



Integrating Eqs. ((3T|) and (|32p tcrm-by-term. we can 
calculate how much energy and angular momentum are 
radiated in each of the dominant modes, similar to the ap- 
proach of Ref. [tHI • We introduce the quantities E^rn and 
Jirn as the total energy and angular momentum radiated 
in each {£, m) mode, computed by integrating Eqs. (pij) 
and (j32p in time, term by term (for conciseness, we com- 
bine both the TO and —to terms into E^rn and Jim and 
restrict our notation to m > 0). Note that while Eim is 
always positive, J^n can also be negative, corresponding 
to angular momentum in the —z direction. These results 
are shown in Table [TTl along with the contributions from 
just the RD phase {t > tpoak, where ipoak is the point 
at which |/^^| reaches its peak, closely corresponding to 
the peak in GW energy emission). We will see below in 
Section |V] that these various energy contributions agree 
closely with the Newtonian predictions for the relative 
mass-scalings. For example, the energy E22 in the inspi- 
ral phase should scale as 77, while the RD contribution 
should scale like 77^. It is important to note that the 
different moments have different scalings: i?33 ~ rfSm?, 
while the /^^ contribution has a much weaker dependence 
on mass ratio: i?44 ^ rf'{l — irj)^ . 



the angular momentum in the inspiral is 

i-U, 



J2 



1 
8^ 



(2)7-22* (3)7-22 



(33) 



=0 



As we will see below in Section |Vl for all the other 
energy and angular momentum modes, the fluxes from 
Eqs. ([31]) ■ ([32|) scale as ui^^^^ or higher powers, and thus 
converge when integrated over uj^^^^^^du). 



IV. MULTIPOLE ANALYSIS OF THE 
NUMERICAL SIMULATIONS 

In this Section we want to investigate how the different 
multipole moments evolve during the inspiral and ring- 
down phases of BH binary mergers. 



A. Inspiral phase 



In the limit of very large initial separation (small ini- 
tial frequency) , each of the Ef„i and Jtm should converge 
to a finite value, with the notable exception of J22. It is 
well-know that the angular momentum of a binary sys- 
tem scales as i?^/^, and is thus unbound in the limit of 
i? — > 00, but it is interesting to see that the higher-order 
contributions to the angular momentum all converge at 
large R. This can be understood directly from Eq. (1321) in 



the Kcplerian limit of i? = AI^^^uj '^1'^ . At leading order, 
radiation reaction follows the relation dt ^ uj^^^^^dcj so 



As can be derived in PN theory [3y| and has been 
confirmed numerically in Refs. [3, Q, the i = 2,m = 2 



mode in Eq. pop is circularly polarized to leading order 
throughout the coalescence. Because of this, Ref. [t^ 
defined the (dominant) orbital angular frequency as 



(34) 



Here, we extend Eq. p4p by defining several (dominant) 
orbital angular frequencies, each of them being related to 
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a specific multipolc moment, J^™ or S^"^, as 



be supported further by the analysis in Sec. IVI Al 



,Ihn 




1 / Q^rn 



(35) 

We plot these frequencies in Fig. [5] for the dominant 
multipole moments P"^, P^, I^"^, and S"^^, for the 
NEgif (left panel) and NEjjf (right panel) runs. The 
amplitudes of the and J^^ modes are too weak and 
dominated by noise to extract a dominant frequency. In 
this figure, as well as most shown in the rest of the pa- 
per, we plot the time variable with respect to ipcak- We 
notice that the frequencies corresponding to the modes 
with i = m agree quite well throughout the inspiral 
and ringdown, but the frequency of the S^^ mode de- 
couples from the others approximately 5QM before the 
peak in the /^^ mode. As we shall see in Sec. I VII this 
is due to the fact that, during the ringdown phase, the 
dominant angular frequency associated to the mode 
is almost twice as large as those of the other leading 
modes (t^, [zH, II^]- This decoupling plays a major role 
in determining the shape of the kick and anti-kick (see 
Sec. I VII below), and also suggests that the transition to 
RD may begin long before the peak of the GW fiux. Sim- 
ilarly, the S^^ mode should converge to a higher RD fre- 
quency ([J320/2 ~ 0.37/Ajft for these runs), but may be 
limited by numerical noise here, as well as possible mode 
mixing with the dominant I^^ moment. 

In Fig. [3] we show the amplitudes of the multipole mo- 
ments in Eq. ([2T|) . Again, the left panel refers to the 
NE§5^ run, while the right panel to the NEj5^ run. The 
mass-quadrupole moment I^^ clearly dominates in both 
cases, while the I^^ and I^^ modes are so weak as to 
be almost completely overwhelmed by numerical noise. 
In addition to having dissimilar amplitudes, the different 
moments also peak at slightly different times, which may 
be related to the fact that RD modes are excited at differ- 
ent times. In particular, the modes mentioned above with 
i m tend to peak later in time, perhaps due to a longer 
transition to the higher QNM frequency. As we shall see 
in Sec.|Vl as the mass ratio becomes more extreme (i.e., 
decreasing 77) , the higher-order modes increase in relative 
amplitude, with I^^ and 5*^^ both proportional to ry 5m. 
I'^^ and S^^, however, scale as 7y(l — 877), so they increase 
only slightly in the range of masses considered here. 

Next, in Fig. [H we show the amplitude of the lin- 
ear momentum flux from the mode-pairs included in 
Eq. (pij) . Here we define the complex flux i^2i,22 _ 
(-14i/3367r)52i/22* and other i^^™.^'™' analogously 
from Eq. (|2ip . As in Fig. O the mass-quadrupole terms 
dominate, with significantly smaller contributions from 
the S^^ and I^^ modes. However, note the appreciable 
flux amphtude from the f33'44 ^ j33j44* ^^^^^ which 
is formally a higher-order correction in a (1/c) expan- 
sion [4^, . From Fig. |4l we expect that the first three 
pairs of modes in Eq. (|2ip should contribute most signif- 
icantly to the recoil. Including the complex phase rela- 
tions between the different modes, we find this result will 



B. Ringdown phase 

We now extract the QNMs, notably the fundamental 
and the first two overtones, present in the most signifi- 
cant multipole moments during the RD phase. We fol- 
low the procedure outlined in Ref. [t^. To avoid possi- 
ble constant offsets introduced by integrating Eqs. (|30ap . 
(|30bp . we prefer to extract the QNMs directly from the 
-2Cim instead of using Iim or Sim- Additionally, from 
Eqs. (I30a|, (l30bl) . we see that (^)/^" and (^)S'^'" are made 
up of both -2Cim and _2C^-m, which in general do not 
have the same QNM frequencies, so it is more reliable to 
extract the RD modes from just _2C^m (however, in prac- 
tice we find that the RD phase is dominated by modes 
with positive m). Following the approach of Ref. [t^, we 
define the complex frequencies uimn'- 



air, 



(36) 



and each RD mode is proportional to exp(— icr£,„„t). In 
this notation, LOimn are the QNM oscillation frequen- 
cies [not to be confused with the dominant frequencies 
of Eq. psp ] and r^mn are the mode decay times, all func- 
tions of the final black hole mass and spin. The sub- 
scripts £ and m are the same spherical wavenumbers used 
above, and n = denotes the fundamental mode, with 
71 = 1, 2, • • • , corresponding to the higher overtones. The 
fundamental QNM frequencies (Jimo are listed in Table Hill 
for the NR runs listed above. All frequencies and decay 
times are measured in units of the final mass Mf . 



We present the RD analysis only for the NEg[f run, 
but the others are qualitatively very similar. We have 
extracted the various QNM contributions to the _2Cfm 
RD signal in the following way (see also Ref. [l^): We 
expect that at late times the 77 = QNM dominates. We 
fit the signal after time tpcak+ir to this single mode using 
non-linear regression and choose tr to minimize the error 
in the fit. We have four dimensionless parameters in this 
non-linear fit: the QNM amplitude and phase, C^mo and 
4>imO, and the QNM frequency and decay time Muimo 
and Timo/M. However, instead of fitting directly for 
these four parameters, we treat Mujgrno and Tgrno/M as 
functions of Of/Mf and Mf/M (which can be obtained via 
interpolation from tabulated values given in Ref. fl%). 
The advantage of using (a{/M{, Mf/M,CimOi(t'emo) for 
the set of fitting parameters comes when we fit to higher 
overtones. As done in Ref. [t^, we extract the QNMs 
treating the real and imaginary parts of -2Cem as in- 
dependent. Below we shall list results obtained from 

Re[_2C£m]. 

By applying this procedure to the dominant mode, 
„2C22, we obtain a{/M{ = 0.669 and M/Mf = 0.965 to- 
gether with the amplitude and phase of the fundamental 
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(t-tpeak)/M (t-tpeak)/M 



FIG. 2: Dominant orbital angular frequency obtained from the individual radiative multipole moments, as determined by 
Eq. (|35}. The different frequencies with £ — m agree closely throughout the inspiral and RD phases. The frequency with 
£ — 2,m = 1 decouples from the others at earlier time and reaches a much higher plateau. The left panel refers to the NEgif 
run and the right panel to the NEjff run. We denote with tpcak the time at which I^^ reaches its maximum. 




(t-tpeak)/M (t-tpeak)/M 

FIG. 3: Amplitudes of the dominant radiative multipole moments. On the left panel we show the modes for the NEoif run, 
while on the right panel the modes for the NEJq^ run. The leading-order mass quadrupole I^^ is about an order of magnitude 
stronger than any other mode. The oscillating behavior of the S^^ moment during RD is likely due to mode mixing with I^^. 
We denote with tpeak the time at which I^^ reaches its maximum. 



TABLE III: Frequencies and decay times for the fundamental QNMs for each of the numerical simulations, ujemo is in units of 
Mr^ and remo is in units of M{. 



Run 


at/Mi 




T21() 


<J.'220 


T220 


<^320 


''"320 


^330 




<^440 


T440 


EQ+_ 


0.697 


0.454 


12.2 


0.531 


12.4 


0.758 


11.9 


0.841 


12.0 


1.14 


11.8 


NEgo' 


0.675 


0.450 


12.1 


0.521 


12.2 


0.749 


11.7 


0.827 


11.9 


1.12 


11.7 


NEjo' 


0.633 


0.442 


11.9 


0.505 


12.1 


0.734 


11.6 


0.803 


11.7 


1.09 


11.5 


NEjb* 


0.423 


0.411 


11.5 


0.445 


11.5 


0.674 


11.1 


0.711 


11.1 


0.963 


10.9 




0.640 


0.443 


11.9 


0.507 


12.1 


0.736 


11.6 


0.806 


11.7 


1.09 


11.5 




0.704 


0.456 


12.2 


0.533 


12.4 


0.760 


11.9 


0.845 


12.1 


1.14 


11.9 
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QNM. We include additional overtones (n > 0) succes- 
sively. For each value of n, we refit the entire function, 
so for n = there are 4 parameters in the fit, for n = 1 
there arc 6, for n = 2 there are 8, and so forth. Thus, ap- 
plying a 6-parameter fit we successfully extract also the 
first overtone simultaneously, obtaining slightly different 
values for Of/Mf = 0.661 and M/Mf = 0.958. Wc find it 
impossible to extract, with a single 8-paramctcr fit, also 
the second overtone. By contrast if wc keep Of/Mf and 
M/Mf fixed and equal to the values obtained when ex- 
tracting the fundamental QNM, we find that we can fit 
up to the second overtone. Moreover, quite interestingly, 
the fit provides waveforms that compare very well with 
the NR waveforms up to the peak of J22 , as can be seen 
in the upper left panel of Fig. [51 

The remaining panels in Fig. [5] show results for the 
other relevant modes -2C33, -2C44 and -2C32- As ob- 
tained in Rcf. [7^, we find a "mode- mixing" in -2C32, 
i.e., the RD waveform is a combination of ^ = 2, m = 2 
and i = 3,m = 2 QNMs. This effect appears to be most 
important between modes with the same m value, and 
may possibly be explained by the fact that the QNMs 
should really be expressed as spheroidal, not spherical 
harmonics [tI, [tI] . Including both sets of modes means 
that the _2C'32 is actually fit using 14 parameters: the 
final mass and spin, and the amplitude and phase of 6 
QNMs. 

By fitting the fundamental QNM for each ringdown 
waveform, we obtain Of/Mf = 0.671 and M/M{ = 0.972; 
af/Mf = 0.527 and M/Mf =^ 0.884; af/Mf = 0.686 
and M/Mf = 0.981, for -2C33, _2C44 and -2C32, re- 
spectively. We also are able to extract the fundamental 
QNM for the -2C21 mode (not shown in Fig.[S]) and find 
af/Mf = 0.678 and M/Mf = 0.960. AU of these values 
for the inferred final BH spin and mass are rather con- 
sistent, except for -2(744. This discrepancy might be due 
to numerical resolution effects, and will be the object of 
future investigations. 

Thus we find that although wc cannot simultaneously 
extract three QNMs (the fundamental and two overtones) 
and wc arc not able to clearly determine the onset of the 
RD phase, we do obtain that for t > fpcak the numer- 
ical waveforms can be well fitted by a superposition of 
three QNMs. This result explains why the simple match- 
ing procedure from inspiral to RD adopted in the EOB 
mod el I39l . l47l .f73t can almost always work succesfully (see 
Ref. [73] for some caveats). In Scc. lVBl wc shall adopt the 
same matching procedure of the EOB model when build- 
ing the full waveform using the pseudo-analytic model of 
Scc.lVl 



V. EFFECTIVE NEWTONIAN MODEL 

In an attempt to better understand the amplitudes and 
frequencies of the various modes during the inspiral and 
merger phases, we present here what we call the "effec- 
tive Newtonian" (cN) model. It begins with calculating 



the leading-order Newtonian formulae for each multipole 
moment of the source, as a function of the BH masses, bi- 
nary separation R, and orbital phase cj). To extend these 
formulae through the end of the inspiral and into the 
merger phase, we introduce an effective radial separation 
to absorb PN effects into the leading-order multipole ex- 
pressions. Each multipole moment is then individually 
matched to a linear superposition of ringdown modes, 
as is done in the effective-one-body model [H, 113, [ll]. 
Taken together with the match to Kerr QNMs, this eN 
model provides an excellent framework within which we 
can understand the details of the linear momentum flux 
and net recoil velocity. 



A. Newtonian Multipole Moments 

Working at leading Newtonian order for each mode, 
we equate the radiative multipole moments to the source 
multipole moments. Restricting ourselves to circular, 
planar orbits, wc find that for non-spinning systems, the 
dominant modes are 
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(37e) 
(37f) 
,(37g) 



where R is the radial separation and co = (j) is the binary 
orbital frequency. Considering only the mass quadrupole 
terms in the linear momentum flux (i.e., the terms pro- 
portional to 5'2i/22*. /3i/22*^ and P^P^*), we obtain the 
well-known result valid at Newtonian order 14711: 



-^105 M^ ^^^ 



7 J4 



(38) 



Including the next-highest order moments in Eq. (jl9[) , we 

get 



.11120(5™ 
' 1323 M 



(39) 



While there may also be next-to-leading order contribu- 
tions from a PN expansion of the multipole moments in- 
cluded in Eq. (fT7)) that would show up in Eq. ([55)1 . we can 
effectively absorb those corrections into the R variable, 
as will be described below. 
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FIG. 5: Comparison of numerical and QNM waveforms for the NEo[f run. The dominant modes analyzed are -iCii {upper 
left), -2C33 {upper right), -2C32 {lower left), and _2C44 {lower right). Note that the -2C32 waveform includes contributions 
from the £ = 2,m = 2 modes as well. We denote with Jpeak the time of the peak of 7^^. 




FIG. 6: Effective radius for different modes, derived from Eqs. (|35|l . (|37ap -( [37g| . The close agreement for the suggests 
we can use a single effective radius Rcs{t) for the Newtonian expressions. We believe that the large oscillations in R^g are 
due to initial eccentricity at early times. Also plotted is the ADM radius (dashed curves) derived from the orbital frequency 
via Eq. (|4ip . the coordinate separation of the BH punctures (dot-dashed curves), and the empirical fit Rat (dotted curves) 
obtained by shifting -Radm by 0.65. The results correspond to the NEoij^ (left panel) and NEoii^ (right panel) runs. We denote 
with tpcak the time at which I^^ reaches its maximum. 
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Combining Eqs. ([38)) and ((39)) wc find the linear mo- 
mentum flux scales like 



Sm 



M 
3 dm 
2 m" 



3475 



l + ^(l-37y)i?2t^2 



1827 
0.977), 



(40) 



which is remarkably similar to the result found in Ref . Q . 
Here we have used R^oj^ w 0.23 — 0.25 at the peak of the 
energy flux, which seems to be quite robust for a range 
of mass ratios. However, the extremely close agreement 
with Ref. is probably to some degree a coincidence, 
since this simple Newtonian formula does not include any 
details of the phase relations between different modes, 
which become especially important during the transition 
from inspiral to ringdown (see Sec. IVI Bl below). Since 
Eq. (|40)) really only applies to the inspiral portion, if 
anything, it should be a predictor of how the peak recoil 
velocity scales. This is not necessarily the same as the 
final recoil, since we find that more extreme-mass-ratio 
BH binaries have a relatively smaller anti-kick, which 
should also play an important role in the scaling relation 
of Ref. @. 

If we compute the above multipole moments p7a[) - 
P7g[ ) using u! as given by Eq. (|35p and R as obtained 
from the puncture trajectories, we do not find a very good 
agreement with the numerical results. This is not sur- 
prising since there is no reason to believe that the New- 
tonian approximation should work well all along the in- 
spiral phase. We should expect that higher-order PN cor- 



rections become important as we approach the merger. 
Furthermore, i? is a coordinate-dependent quantity, and 
thus does not necessarily have the same meaning in a 
PN expression as in NR. Since our scope is limited to a 
diagnostic of the NR results, and not to a precise com- 
parison with PN calculations, instead of including PN 
corrections in Eqs. (|37| - (|39)) . we investigate whether by 
properly scaling the Newtonian expressions we can get a 
better agreement until the merger. We can also think of 
this normalization as a way of resumming the PN expan- 
sion. 

Quite interestingly, if we compute the amplitudes |/^™| 
or 1 5*^™ I from the numerical data, and the angular fre- 
quency uj from Eq. (|35)) , we find that the radii i?^™ which 
appear in the RHS of Eqs. (|37ap -( [57g| ) are rather inde- 
pendent of the multipole moments £ and m, as Fig. [6| 
shows. We denote the radii i?^™ computed numerically 
as effective radii i?^ff ■ "^'^^ close agreement between the 
frequencies (see Fig. [2|) and effective radii for each mode 
suggests we can use the Newtonian expressions and a 
single i?off(i) and orbital frequency (jj{t), e.g., R^{t) and 
Lo^^iox all modes with a high degree of accuracy for the 
entire inspiral phase and even during the transition to 
merger. 

For comparison we also show in Fig. [5| the radius 
from the puncture trajectory (dot-dashed curves) and 
the radius computed using the Arnowitt-Deser-Misner 
transverse-traceless gauge (dashed curves), given as a 
function of frequency through 3PN order by [74| 



i?ADM = Af'/'c^"'/' 



l+UJ 



^'=(--I)-'"H-?"4)-'( 
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77 TT 



81 



(41) 



Here we use the orbital frequency oj derived from the 
I^^ mode via Eqn. ()35p . giving a constant value during 
the RD phase when the orbital frequency is meaningless. 
Fig. [6| shows interesting agreement between i?ADM and 
the radius from the puncture trajectory, and a constant 
offset between i?ADM and i?eff- The latter is due to the 
fact that the amplitude of the multipole moments com- 
puted at leading Newtonian order does not reproduce the 
numerical relativity amplitude [gl, [t^ , and higher order 
PN corrections need to be included. Motivated by this 
similarity between i?ADM and i?cff , we attempt to fit em- 
pirically the i?eff curves in Fig.[B|by simply shifting i?ADM 
by 0.65. The fit curve is included as a dotted curve in 
Fig. [6l As we accumulate longer and more accurate NR 
data for a wider range of 77 values, and study possible an- 
alytic resummation of higher-order PN amplitude correc- 
tions, we should be able to work out a widely applicable 
amplitude-scaling factor to be included in leading-order 



analytic waveforms [77|. 

In the next section, we shall investigate how this sim- 
ple eN model can be combined with a superposition of 
QNMs, as described in Sec. HVBl giving a good repre- 
sentation of the NR results. 



B. Matching to ringdown 

We now match the inspiral and RD waveforms in a 
mode-by-mode fashion following the philosophy of the 
FOB approach [3§| . Note this is not the same analysis of 
Section llVB) where we fit the numerical data throughout 
the RD phase with a superposition of QNMs. Here we 
match the data at a single point at the transition from in- 
spiral to RD and see how well it agrees with the rest of the 
RD phase. A similar attempt was followed in Ref. [IJ, 
where for simplicity the authors performed the match- 
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ing to the Schwarzschild QNM frequencies, while we use 
the Kerr QNM frequencies and match to the fundamen- 
tal QNM frequency and the first two overtones, as done 
in Ref. [zl]. We olatain the QNM frequencies and decay 
times from Ref. [zH as a function of af/A'/f (taken from 
Table H] above). For the fundamental and two overtone 
QNMs, we can match a given multipole mode by equat- 
ing it and two time derivatives to a linear combination 
of QNMs. 
We write 



/^"(i):=A(t)e-'^(*) =^A,„„e- 



■iO'f„,„(t — t„atch) 



(42) 



where the complex QNM frequencies are known functions 
of the final BH mass and spin, and we must solve for the 
complex amplitudes Ai,nn ■ Matching three QNMs we get 



I (^match) — ^ ^ Almn, 



(43a) 



— /^™(<„iatch) — —i (^Imn^emn, (43b) 



n=0 
2 



^/^'"(tmatch) = -y^q-^mn Ajmn, (43c) 



or as a simple matrix equation 



-G 
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1 jlm 








Atm2 1 







(44) 

In Fig. \7\ we compare the NR modes to the modes 
obtained by the effective Newtonian model described in 
Sec. IV Al until imatch and by the superposition of three 
QNMs for t > tmatch- During the inspiral, the different 
moments are calculated according to Eqs. (I37ap -( [57gl ), 
using a single i?off and lud determined from the I^^ mode, 
with the exception of the mode, where we instead 
use the higher frequency w^^^ (but same i?cff)- We treat 
^match as SL frcc parameter: if we stop the inspiral too 
early, the eN mode amplitudes are still growing, so the 
sudden transition to decaying RD modes prematurely re- 
duces them. On the other hand, if the inspiral is contin- 
ued too long, we tend to lose the important phase shifts 
between the modes that only begin during the transi- 
tion to RD. This is particularly evident in the I'^^ mode, 
which undergoes an unexplained phase-shift around the 
transition to RD, and also decays at somewhat differ- 
ent rate than is predicted from QNM theory (see above, 
Sec. lIVB|) . Motivated by the results of Sec. lIV Bl notably 
by the fact that a superposition of three QNMs can fit 
very well the NR waveforms starting from the peak of the 
energy flux, we choose as best-matching point the peak 
of the energy flux. 

Having shown a reasonably close match for each of 
the radiative multipoles between the effective Newtonian 



model and the numerical data, it stands to reason that 
the total recoil calculated with this model should agree 
as well. This is shown in Fig. [51 where we have also 
varied the matching point around ipoak- We first note 
the close agreement between the eN models with vary- 
ing imatch, suggesting the inspiral-to-ringdown matching 
method described above is relatively robust. Not sur- 
prisingly, since the individual modes agree, we also find 
reasonable agreement between the NR data and the eN 
predictions for the recoil. 

However, this agreement may be partially fortuitous, 
since the eN model cannot predict the mode phase shifts 
around t = tpoak, most notably that of the /^^ mode 
described above. In Section fVl Bl below, we will examine 
this phasing in greater detail and show how it affects 
the overall kick. At this point, we unfortunately do not 
have a clear understanding of the underlying cause of 
the phase shift, but it may well be related to the slightly 
different times of transition from inspiral to ringdown for 
the different modes. Preliminary results also suggest that 
this de-phasing effect is reduced in more extreme-mass- 
ratio systems, as we shall see in Appendix El 



VI. ANATOMY OF THE KICK 

In the above Sections, we have laid the groundwork 
for a multipolar analysis of the gravitational recoil, de- 
scribing the momentum flux as a combination of radiative 
multipole modes. Along with the psuedo-analytic models 
for the inspiral and ringdown phases, we can now give a 
detailed description of the "anatomy" of the kick, namely 
the way the different modes combine to produce a peak 
recoil velocity, followed by a characteristic anti-kick and 
then asymptotic approach to the final value of the BH 
recoil. 



Contribution from different moments 



In Sec. nil Al we showed how the radiative multi- 
pole moments contribute to the linear momentum flux 
through the integral of the *4 scalar [Eqs. p0)) . p2)) ]. 
Here, we want to determine exactly which modes we need 
to include in the multipole expansion Eq. to get a 
good representation of the full recoil, and which are the 
pairs of modes in Eq. pT|) that contribute most. 

By including only a select choice of terms in the ■04 ex- 
pansion Eq. pop , we can calculate the linear momentum 
flux by direct integration of Eq. (|12[) and compare it with 
the predictions of Eqs. (fT7)) - pT|) . in each case including 
only the appropriate moments. This is a good way of 
double-checking those lengthy equations term-by-term, 
and in practice we find excellent agreement, limited only 
by the numerical accuracy of the simulations. Similarly, 
we can use this method of truncated expansion to deter- 
mine which modes are necessary for calculating the recoil 
up to a given accuracy. The results of using higher and 
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FIG. 7: Comparison of the effective Newtonian and NR radiative modes during inspiral, merger and RD phases. The data 
refer to the NEo5^ run. We denote with tpoak the time at which I^^ reaches its maximum. 



higher order multipolar moments are shown in Figs. [H] 
and [TOl for the NEgj]^ and NEgj,^ runs, respectively. 

In the left panels of Figs. [9] and [10] we show with a 
solid curve the exact recoil velocity from Eq. (fT2|) . with 
a dashed curve the contribution from terms up to £ = 4, 
i.e., those obtained from Eq. (fT7|) and and with a 

dotted curve the contribution from just the three lead- 
ing terms in Eq. (HI]), valid for non-precessing BHs with 
kicks in the orbital plane. We conclude that the linear 
momentum flux is dominated by the p^p"^*^ j33j44*^ 
and 521/22* tej-jns, which combine to produce the pri- 
mary kick and anti-kick agreeing with the exact result 
within < 10% throughout the entire merger. Note that 
the flux from the S^^P^* term, while not insignificant 
in Fig. m contributes almost nothing to the net recoil 
velocity. This is largely due to phase relations between 
the various modes during the transition from inspiral to 
ringdown, described below in Sec. IVI Bl 

In the right panels of Figs. [S] and [TO] we show the dif- 
ference between the calculation obtained including terms 
up to £ = 3,4,5,6, and the exact result. It seems 
clear that we need modes up to and including £ — A 
to get an accurate estimate of the recoil velocity. For 



more extreme mass ratios, higher-order moments be- 
come relatively more important, but remain strongly sub- 
dominant to the £ < 4 modes [ill . [TSj . 



To understand more clearly the relative contributions 
of the different modes to the total recoil, we will in- 
clude analysis of a few more simulations including non- 
precessing spins. As mentioned above in Sec. IIII Al non- 
precessing spins do not introduce any additional mo- 
ments compared to the non-spinning simulations, but 
simply modify the relative amplitudes of the different 
modes in Eq. (|2ip by adding the spin terms. Thus, once 
we determine how the spins modify the individual modes, 
we can use the same analysis for the spinning and non- 
spinning cases. 



Again equating the radiative multipole moments with 
the source moments, we get the leading order spin-orbit 
modifications to Eqs. ([T7a|)-([57g]) [see Eqs. (3. 14), (3. 20) 
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where we have introduced the spin vectors 
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FIG. 8: Comparison of the effective Newtonian model and NR 
predictions for the recoil velocity for a range of inspiral-RD 
matching points. We denote with fpoak the time at which I^^ 
reaches its maximum. The data refer to the NEoif run. 



in Ref. and Eq. (5.5) in Ref. [Zi]]: 
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In all of the simulations considered here, the dimen- 
sionless spins are equal (|ai|/mi = |a2|/m2) and point 
in opposite directions, = 0, so for the leading-order 
terms in Eqn. (|2ip we are left only with the modifica- 
tions of S*^^ and I^^, due to and S33, respectively. 
Then Eqs. (|37|1 and (|45p give the linear momentum flux 
during the inspiral for each of the first three dominant 
terms in Eq. 0J^ : 



= ^^^R'u;H2SmR'u. + 3A^)e^t (47a) 
F'lt' = -y^f^(l-3r/)i?^c.^5m-fc.E^3)e'*. (47c) 



While these flux formulae contain terms of various orders 
in w, we expect that the effective Newtonian scaling of R 
ensures that we are including all relevant PN terms, at 
least in the cases where the Sm terms dominate over the 
spin corrections. When the spin terms begin to dominate, 
we find that it becomes more difficult to use a single 
effective R for all modes. This can be seen in Fig. [TTl 
which plots i?cff as in Fig. [6l but for the NE^i^ run, 
where the and Sm terms in Eq. (j47ap arc comparable, 



making it difficult to derive a reasonable i?efi(S'^^). 

Even for non-spinning runs, in order to get reasonable 
agreement with the NR data, we find that one must be 
careful towards the end of the inspiral to distinguish be- 
tween u^'^ and u^^^ in Eq. (|47ap : 

^l^sf « (MVM)i?^ (c^^^)^ {u;^'^ {2SmR'u;^'^+3An . 

(48) 

The amplitudes of these fluxes are plotted in Fig. [T^] 
for the four runs NE2:3_, NE2;3 , NE^jf^ ^nd EQ+_. As 
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FIG. 9: In the left panel we show the net recoil kick, integrated from the linear momentum flux via Eq. (|12|l (solid curve), from 
all modes with £ < 4 (dashed curve) and also limiting the modal composition of ^4 to just the three dominant mode pairs in 
Eq. (|21|l (dotted curve). In the right panel we show the difference between the exact result and the ^4 expansion Eq. pUjl . 
limited to ^ < 3, 4, 5, 6. The data refer to the NEnrf run. We denote with tpoak the time at which I^^ reaches its maximum. 
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FIG. 10: Same as Fig. U but for the NE^ run. 



seen in Table U the NE^:^ run has = 0.2A'P, while 
the NE^-l run has = — 0.2M^, respectively adding 
destructively and constructively with the Sm term in 
Eq. (|47a|) . This difference is clearly seen in the blue 
curves in the top two panels of Fig. [T^l Also no- 
table in these plots is the somewhat smaller difference 
in the amplitudes of ^22, 33^ ^^^g similar effect from 
the constructive/destructive additions of Sm and SI3 in 
Eq. I|47bp . As we see in Fig. [TH NE^ff appears to be 
the average of NE^'i and NE^'^, while the flux from 

EQ^ is strongly suppressed due to the Sm = terms 

in Eq. (j47p , leaving only the flux from the terms propor- 
tional to A^ = -0.2M2 and Efg = 0.075. However, as 
noted above, when the spin terms dominate the flux, as 
in the case of equal-mass BHs, the eN model with a sin- 
gle i?eff begins to break down. Yet even in this situation, 
Eqs. (|47ap - (|47c[) still have qualitative (if not quantitative) 
predictive value, including the relative phases between 



the different mode-pair fluxes during the inspiral. 

In each panel of Fig.[T21 we also plot with dashed lines 
the eN prediction for the various flux amplitudes. In 
almost all cases, the eN flux is quite close to the NR 
results up to about lOM before ipoak, when the eN model 
begins to break down, especially for the spinning runs. 
The amplitude differences near the peaks are comparable 
to those seen in Fig. [7] for the NEpij^ run. The notable 

exception is the F^^'^^ flux from the NE^i^ and EQ^ 

runs, where the spin terms dominate over the Sm terms. 



B. Transition to ringdown and the de-phasing of 
the multipole modes 

Since the flux vectors defined by Eq. (|47p will not gen- 
erally be co-linear, to understand the time evolution of 
the recoil velocity, we must first understand the phase 
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FIG. 11: Rcfi derived from different multipole modes, as in 
Fig. [61 for the NE!;+ run. The 5^^ mode for this run has 
comparable contributions from 5m and A', making it difficult 
to derive a reasonable Rcs{S'^^). 

relations between the different modes. From Eqs. ([37]), 
and (|47p . wc see that during the inspiral phase, the 
individual moments and the resulting flux vectors evolve 
according to a single orbital phase 4>, with F^^^^^ pointing 

in the opposite direction to F^^^^^ and F^^^'^^. However, 
as we can see from Fig. ^ as the binary evolves from in- 
spiral to RD, the frequency (and thus phase) of the 
mode decouples from the other dominant modes. Upon 



closer inspection, we find that even the /22 /33 g^j^^ j44 
modes deviate from each other enough to undergo a sig- 
nificant phase shift at the inspiral-RD transition. 

To quantify these effects, we define the following phase 
differences: 
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Here we use the notation ^/i™"™ to describe the phase 
difference between two complex flux vectors, where m 
and to' correspond to the larger w-values of each mode 
pair that makes up the flux. These deflnitions are valid 
throughout the inspiral, merger, and ringdown phases. 
In the inspiral phase, we can see that for the unequal- 
mass runs where Sm dominates with respect to the spin 
terms in Eqs. (|47ap - (l47cp . we have 

cos ^Lp = COS = - 1 , COS V-f-p" = 1 . (50) 

For the EQ+_ run with 5m ~ 0, Eq. ([T7)) predicts that all 
phases have cos V'insp = 1 during the inspiral (as shown 
in TablelH and Sfg have opposite signs, so all the flux 
vectors in Eq. (|47)) are parallel). During the RD phase, 
using Eq. ([32), we can approximate the flux vectors and 
phase evolution in terms of the fundamental QNM fre- 
quencies aira^: 



p21.22 

^RD 

p22,33 

^RD 

p33,44 
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^matcli exp[-i(cr210 - (y*21^){t - ^match)] 
i(o'220 ^ cr33o)(^ ^ ^match)] 
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where the F^^^^ fluxes include complex phase information at the matching point. Taking the phase differences 



between these RD modes gives 

cosiArd^ 
cosV'rd^ 
cos^'Id^ 



COs[(lJ210 - 2^220 + t^330)(i ~ tmatcli) + ^m^Llil ' 
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COS[(W220 - 2^330 -|- W44o)(t - tmatcli) + ^f„'^tclJ ' 



(52a) 
(52b) 
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Here <I>match is a phase offset determined at the tran- 
sition from inspiral to ringdown. Quite interestingly, 
we find that for the range of final BH spin parameters 
0.5 < Of/A/f < 0.8, the linear combinations of frequen- 
cies in Eqs. (|52ap - (|52cp vary by less than ~ 5%. Thus, 
if we compute the above expressions for the a;/,„o corre- 



sponding to flf/Aff = 0.7, we have [72 
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FIG. 12: Relative amplitudes of the dominant multipole mode-pairs in the linear momentum flux. Also shown in the dashed 
curves are the eN model predictions for the flux amplitudes. We denote with tpeak the time at which /^^ reaches its maximum. 



Even more intriguing, we find that for the unequal- 
mass simulations described above, the phase relations 
during the inspiral and RD are almost identical, regard- 
less of spin orientations. This can be seen clearly in 
Fig. 1131 which plots cos')/' during inspiral, merger and 
RD for the different runs. The colinearity of the flux 
vectors is clear during the inspiral phase, and the sinu- 
soidal oscillations of the phases during RD agree well 
with the analytic predictions (plotted in dashed curves 
in Fig. [T3)) . Since the analytic models are most reliable 
during the inspiral and RD phases (but have more dif- 
ficulty tracing the merger portion), we omit in Fig. [T3] 
the transition region of —IDA/ < (i — ipoak) < lOM. 
The analytic phase relations during inspiral are deter- 
mined by Eq. (jSP]) and during ringdown by Eqs. (|53ap - 
(|53c[) . Here we use a tmatch (and corresponding <i>match) 
about 20M after fpcak to ensure that the multipole mo- 
ments arc truly dominated by the fundamental QNMs, 
and thus Eqs. (|53a[) - (|53cp are valid. Note that the phase 
differences for EQ-| arc particularly noisy since the am- 
plitude of the l^^ moment is zero to leading order, and 
thus it is more difficult to extract a clear phase for that 
mode. 



The feature that is most difficult to explain from an 
analytic model alone (and is thus omitted from the eN 
curves in Fig. ll3p is the roughly 180-degree jump in phase 
between F^^^ and i^j^^p^, beginning around 20 A/ before 
the peak. This appears to be a feature in all the unequal- 
mass runs examined, but preliminary results suggest that 
is less significant (i.e., a smaller phase shift) for more 
extreme-mass ratio systems, as we shall discuss in Ap- 
pendix [XI We arc not able to explain it with the ad- 
ditional RD overtone modes described in Sec. IV Bl but 
using slightly different RD matching points for the dif- 
ferent multipoles may help explain the issue. 

C. The anti-kick 

These fiux amplitudes and phase relations can now be 
used to understand the amplitude of the kick and anti- 
kick, by which we mean the difference between the peak 
and the final recoil velocities (see Fig. [1] for an example). 
Throughout the inspiral phase, the amplitude and rota- 
tional frequency of the fiux vectors in Eq. (|47p are mono- 
tonically increasing, giving the familiar outward-spiraling 
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FIG. 13: Phase differences between diflFerent mode-pair flux vectors, as defined by Eqs. (|49a|) - (|49c|l . The data refer to the 
NE+i (upper left panel), NE?;+ (upper right panel), NEqo^ (lower left panel), and EQ+_ (lower right panel) runs. The dashed 
curves are the eN model predictions of Eqs. (|50|I . (I53|I . We denote with tpcak the time at which I^^ reaches its maximum. 



trajectory for the velocity vector. Then, in the RD phase, 
the dominant frequencies are nearly constant while the 
amplitudes decay exponentially for each mode, giving an 
inward-spiral that decays like a damped harmonic oscil- 
lator around the final asymptotic recoil velocity. 

These trajectories in velocity space can be seen in 
Fig. [T31 along with the instantaneous flux vectors from 
the competing mode-pairs. Clearly, even small changes 
in the mass ratios and spins orientations of the BHs can 
give a rather diverse selection of velocity trajectories. 
Note in particular the difference between the NE^'i^ run, 
dominated by the i^22,33 large anti-kick, and 

the EQ^ run, which in contrast is dominated by the 

^21,22 -y^g gj^^ ^Yia,t the EQ+_ run has no anti- 

kick, which can be explained by the slowly rotating flux 
vector that does not spiral back inwards, but rather drifts 
off slowly towards infinity during the ringdown. The dif- 
ference between these two runs can be explained entirely 



by examining the real part of their fundamental QNM 
frequencies CT^mO: which in turn determine the rotation 

rates of the flux vectors in Eq. ([51]) : EQ-| is dominated 

by CJ220 ~ = 0.08/Mf, a much slower frequency than 
^^330 — 1-1^220 ~ 0.31/A/f, which causes the rapid inward- 
spiral of the NE^ij_ run. 

To calculate the recoil velocity, we must integrate the 
linear momentum flux vectors in time. (For the initial 
velocity vector, we integrate the post-Newtonian approx- 
imation for the momentum flux from t = — oo to the 
beginning of the numerical simulation [46| . This effec- 
tively sets the centers of the spiral curves in Fig. [14] to 
correspond to the origin in velocity space.) We can get a 
reasonable analytic approximation by using Eqs. (|47p and 
([5T|) for the inspiral and RD phases, respectively. In the 
adiabatic inspiral, the complex recoil velocity v = Vx+iuy 
can be written as 
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FIG. 14: The recoil velocity vector evolving in the v^-Vy plane {black solid curve), along with the flux vectors due to the three 
mode pairs at each time interval along the velocity trajectory (same color scheme as Fig. I12|) . The data refer to the NE+l 
(upper left panel), NE?;^- [upper right panel), NEojf {lower left panel), and EQ+_ {lower right panel) runs. We denote with the 
label peak the time at which I^'^ reaches its maximum. 



while for the RD portion we have 



F{t') dt' ~ — Fn,.tcb, (54) 

— oo ^^match 



mD(i)= / F{t')dt' 

' ^match 



^ match 



Vi = WRD(i oo) ~ ^ 



F. 



£ni.£'m 



1^, ^^raa - l^tm'O 



match ' 



(55a) 
(55b) 



summing the contributions from each pair of modes {im, I'm'). Then the total velocity in each of the dominant mode 
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pairs is given by 



^2 ^ 5 2 

^ ^ ^matcli'<^matcli(2'5TOi?„jjjtj,[ja;i„atch + 3A^) 



21,22 



O'210 



'220 



-^niatch '^match('^'^ + ^match^h) 



36 

-y ^ (1 - 3?/) ^match '^match ((^^ + WmatchS^g 

I 



.,22.33 

^t^match e'^'m-tch 



^220 — 0". 
1 



330 



"■330 



'440 



(56a) 
(56b) 
(56c) 



The phase (t^match defined as the angle made between 
the flux vector _F2i,22 ^^^^ ^j^g velocity vector v at the be- 
ginning of the ringdown (with other phases (^^^tch' '^match 
defined analogously). Because of the anomalous phase 
shifts and departure from adiabaticity at the transition 
from inspiral to ringdown, these angles are difficult to 
predict with an independent analytic model, but can be 
calculated easily from plots like Fig. [T31 However, the 
accuracy of Eq. (|56|1 is limited both by the adiabaticity 
condition of Eq. ((54)) as well as the accuracy of the spin- 
orbit corrections to the eN model (see Fig. [TT|) . There- 
fore, in analyzing the anti-kick in terms of RD modes, we 
find it more useful simply to integrate Eq. (j54p directly 
from the numerical data during the inspiral, and then 
attach the fundamental QNM terms from Eq. ([55]) at the 
matching point tmatch = tpoak- 

Given Wmatch at the end of the inspiral, we can use 
this quasi-analytic approach to predict the maximum and 
final recoil velocities (wmax and V{, respectively). These 
predictions are plotted as black dashed curves in Fig. [151 
to be compared with the solid black curves of the exact 
NR results. Within this context, we define the anti-kick 
magnitude as 



ak 



V{ - v„ 



and the net ringdown contribution as 



/] 



RD 



Vf — Wmatch 
^match 



(57) 



(58) 



where Wmax and V{ are the (real-valued) velocity magni- 
tudes calculated analytically from Eq. ([SS]) . 

In the case of the NE^i]_ run, where the recoil is almost 
entirely dominated by the F^^'^^ flux, we find a large 
anti-kick with /ak = —0.53 and /rd = —0.5. On the 
other hand, for the NE^"^ run, as can be seen in Fig. [T51 
the net recoil velocity continues to increase after imatch = 
tpoak before turning around for a small anti-kick of /ak ~ 
—0.11. The total effect of the ringdown phase is actually 
to increase the recoil with /rd = 0.68. An intermediate 
effect is seen for the NEq5^ run, with /ak = —0.28 and 
/rd = —0.04. However, as seen above in Fig. [T31 for the 

EQ^ run, we see no anti-kick, with /ak = —0.01 and 

/rd = 0.58. 



In general, we find the magnitude of the anti-kick is 
primarily dependent on the relative magnitudes of the 
S'^^ and I^^ moments. When S^^ dominates (e.g., when 
6m and add constructively), the ringdown rotation is 
slow and there is a small anti-kick, whereas a dominant 
j33 ]-QO(jg (e.g., large Sm or no spins) gives a rapidly ro- 
tating ringdown flux and thus a large anti-kick. Further- 
more, from Eq. (|47p . we see that for non-spinning BHs, 
both the S*^^ and I^^ modes share the same mass and fre- 
quency scaling, so the relative size of the anti-kick should 
be roughly independent of mass ratio (see Appendix |^ 
for a caveat). 

We would like a more quantitative picture of how these 
flux vectors add constructively and destructively to give 
the total recoil velocity to support the analytic estimates 
presented above. Using v = J Fdt, we can write 



dt' 



dt 



(v-v) 



(59) 



where v-v = 1. Breaking F up into the contributions 
of the dominant modes as above, and then integrating in 
time gives 



V 



v-F^^^^^dt, 
22.33 ^ U.F^^^^Ut, 



which add linearly to give to total recoil velocity: 



,,22,33 



,,33,44 



(60a) 
(60b) 
(60c) 

(61) 



Note that with these definitions, the i/™'^ ™ are all real, 
but can be positive or negative. These different veloci- 
ties are plotted in Fig. 1151 with the same color scheme 
as in Figs. fT2l and [Til along with the total recoil velocity 
in solid black curves. Also shown in Fig. [T5l is the ve- 
locity 7732^33 (dashed blue curves), defined analogously to 
Eq. ([SOal for the S'32/33* terms. The small contri- 
bution from this mode pair further justifies our focus on 
the more dominant pairs of Eq. ([HJ) and Fig. [H 

In the NE^i|. run, where the modal analysis shows the 
^21,22 g^j^^ ^33,44 terms canceling out, we see that 
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FIG. 15: Relative contributions to the total recoil velocity from the different multipole mode-pairs. J^^/sa* ^^^^ curve) is the 
dominant mode for unequal-mass binary systems, while 521/22* ^j^j^g gm-ye) dominates for spinning, equal-mass binary systems. 
Also plotted is the contribution from the 5'32j33» terms (blue dashed curve), demonstrating its very small contribution to 
the total recoil velocity. For t > tmatch ~ tpoak, we include the quasi-analytic prediction for VB.o{t) (black dashed curves), based 
on the fundamental RD modes from Eq. (|55|) . The data refer to the NE^?. (upper left panel), NE?;i]_ (upper right panel), NEpi)^ 
(lower left panel), and EQ+_ (lower right panel) runs. We denote with tpoak the time at which I^^ reaches its maximum. 



the total recoil velocity (black curves in Fig. fTS)) is almost 
entirely dominated by the i^22,33 g^j, curves). On 
the other hand, for the NE^^ run, the F^^^^^ flux is much 
stronger, adding destructively with the ^^2,33 g^j, (j^j-jng 
the RD. This has the effect of both increasing the peak 
velocity and also decreasing the relative strength of the 
anti-kick, due to the slow rotation frequency of the F'^-^''^'^ 
flux during ringdown, as described above. As expected, 
the NEpji^ run displays behavior intermediate between 

these two extremes. The EQ.) run, however, is entirely 

dominated by the i^2i,22 g^^^.^ ^^^^ ^j^^^g experiences no 
anti-kick, but rather drifts off slowly in a nearly constant 
direction, as seen in the bottom-right panel of Fig. [Ml 



D. Application to non-planar kicks 

One of the most remarkable results from the recent 
renaissance in numerical relativity was the prediction of 



extremely large kicks from equal-mass BHs with spins 
pointing opposite to each other and normal to the or- 
bital angular momentum, producing a recoil out of the 
orbital plane [l^, [ij, [l^, . While this configuration 
can produce recoils of nearly 4000 km/sec, the analogous 

non-prccessing configuation (EQ^ run in this paper) 

gives a kick of only ^ 500 km/sec in the case of maximal 
spin [l3, [m, [13 • The multipole analysis tools developed 
above can be used for understanding and explaining this 
remarkable difference. 

First, we should note that leading-order FN estimates 
of the linear momentum flux during inspiral suggest that 
the discrepancy should be less than a factor of two. For 
example, Eq. (3.31b) of Kidder [i^ gives the spin-orbit 
contribution to the momentum flux for circular, Keple- 
rian orbits as 

16 ^ 

Fso = YgM'M^ifix A + (nxv)(v.A)], (62) 
with n and v being the normalized separation and veloc- 
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ity vectors, respectively. For spins parallel to the orbital 
angular momentum, the term in square brackets has mag- 
nitude A^, while for planar spins, it is 2A^' sin</)A, where 
is the angle between A and n, and A^ is the magni- 
tude of A in the orbital plane. 

Not surprisingly, we get the exact same results from 
the multipole analysis of Eqs. p7|) . P5|) . and one 

new multipole moment: 

Sil = i^^r^RLu^e-^*{A^-^Ay), (63) 
I 



where is the orbital phase of the binary. So in both 
paradigms, we see that, when maximizing over sine/) a, 
the planar-spin orientation should result in a recoil twice 
as large as the parallel-spin case, leaving a factor of 
roughly 4 difference unexplained. 

From Eqs. ((64)) . ([65]) we see that the only relevant 
modes involved should be 5*^^, and 5*^^ (for these 
equal-mass systems the momentum flux is dominated by 
a single mode pair, responsible for > 95% of the final re- 
coil value). In the left panel of Fig.[Tn]we plot the ampli- 
tude of from the EQ.j simulation, along with that 

of a planar-spin simulation EQpianar- AH other binary 
parameters and the initial conditions are the same. Re- 
markably, the mass-quadrupole moments I^^ are nearly 
identical (and dominant) in both runs, and this suggests 
that the energy and angular momentum fluxes are the 
same [see Eqs. (|3T|) . ([32|) ]. This is in fact quite reasonable 
since the total spin of the system is zero in both cases. 
However, we see in the right-hand panel of Fig. (TH] that 
the peak amplitude of the S^^ mode is a factor of ^ 2.5 
greater than that of the S'^'^ mode from the EQpianar and 
EQ^ nms, respectively. 

Yet Eqs. (|15|) . (|63p suggest that these two modes should 
have exactly the same magnitudes, at least during the 
inspiral phase, and presumably during the RD as well, 
since the RD amplitudes are completely determined by 
the mode amplitudes at the matching point. It appears 
from Fig. [12] that S^"^ and S^^ do in fact have the same 
amplitude at early times, but the relatively noisy data 
and short duration of the simulations make it impossible 
to say for certain. If this is the case, one possible explana- 
tion for the sudden remarkable increase in the amplitude 



while on the other hand, the mode is zero for the 
planar spin configuration. Combining these equations, 
we get 



(64) 



(65) 



of S*^^ might be mode-coupling with the dominant /^^ 
mode, as the inspiral phase begins to transition to the 
RD phase. This coupling is analogous to that of S^'^ and 
I^^ described above in Sec. IIVBI an effect that is ap- 
parently only important between modes with the same 
TO- number [73. [73j . We hope to address this question in 
the future with longer simulations to confirm the agree- 
ment at early times, as well as other spin configurations 
that should enhance specific multipole modes and may 
help identify other similar cases of mode amplification. 

Lastly, from the ringdown contribution to the velocity 
[Eqs. ([55]) . ([ro)) ]. we can understand another difference 
between the planar- and parallel-spin orientations. In- 
stead of having two different RD frequencies CT210 and 
(T220 combine to give a slowly rotating flux vector, for 
the planar-spin case, we have two identical RD frequen- 
cies for I^^ and S^^ in Eq. ([65]) . giving precisely zero 
rotation to the RD flux. Furthermore, as the spin vec- 
tor A is precessing faster and faster in a positive di- 
rection around the orbital angular momentum vector, 
even during the inspiral the two modes and 5*^^ be- 
come nearly locked in phase, producing a relatively long- 
duration burst of linear momentum flux in a single direc- 
tion during the merger phase. Combined, these effects 
essentially straighten out the spiral curve in the lower- 
right panel of Fig. [TH providing another factor of ~ 1.6 
of increased recoil velocity for planar spins. 

In Fig. [T7] we show the combination of the above ef- 
fects. In the left panel, we plot the linear momentum flux 
from Eqs. ([M)) . ([S5|) . showing the factor of two increase 
predicted by the Kidder formula and our Eqs. p7|) . (|18|) . 
along with the factor of 2.5 increase in the amplitude of 



and using Eq. p8)) we obtain 



F.^^[-28W''S''n] = |^i?^c.«(A^^sin0-A^cos< 

= ^^i?3^6^Psin0A, 

15 M 

I 
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FIG. 16: left panel: Comparison of the multipole amplitudes I^^ for the two different equal-mass simulations: EQpianar (solid 
line) and EQ+_ (dashed line), right panel: The 5*^^ amplitude from the planar-spins run (EQpianar, solid line) and the 
amplitude from the parallel-spins run (EQ+_, dashed line). We denote with tpeak the time at which I^^ reaches its maximum. 




FIG. 17: left panel: Comparison of the linear momentum flux for the two different equal-mass simulations: EQpianar (solid 
line) and EQ-i — (dashed line), right panel: The total recoil velocity from the planar-spins run (EQpianar, solid line) and the 
parallel-spins run (EQ-| , dashed line). We denote with ipeak the time at which I^^ reaches its maximum. 



S^'^ relative to 5^^. In the right panel, we plot the recoil 
velocity for both runs, which includes the effect of flux 
rotation during the merger and inspiral phases, account- 
ing for another factor of ^ 1.6, giving a total discrepancy 
of z;(EQpi,„,,)/«(EQ+_) « 2.5 x 2 x 1.6 = 8. 



VII. DISCUSSION 

In this paper we analysed several numerical simula- 
tions of binary BH coalescence, focusing on the physics 
of the recoil. We developed tools, based on the multipo- 
lar expansion [s^, HH, [H, [H, [11] , that can be used as a 
diagnostic of the numerical results, and understand how 
the recoil velocity evolves during the inspiral, merger, 
and ringdown phases of the coalescence. 

We wrote explicit expressions for the linear momentum 



flux expressed in terms of radiative multipole moments 
through £ = A, valid for generic spinning, preccssing BH 
binary systems. We found that these formulae are suf- 
ficient to obtain the total recoil velocity with high ac- 
curacy. By comparing the amplitudes of the different 
multipole moments, we found that in the case of non- 
precessing spins-and thus a recoil in the orbital plane- 
only three pairs of modes contribute to most of the linear 
momentum flux, notably S^^P'^*, /22/33* ^nd /33/44*^ 
Those modes account for the total recoil with an accuracy 
on the order of ~ 5 — 10% throughout the simulations, 
(see Figs.mni). 

The way in which the contribution from these three 
pairs of modes builds up is not trivial, since not only the 
relative amplitudes, but especially, the relative phases are 
also quite important. We found that the relative phases 
between the three mode-pairs are nearly constant during 
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the inspiral phase, but start diverging at the onset of the 
transition from inspiral to RD (see Fig. fO)) . The late- 
time evolution can be described reasonably well with ana- 
lytic formula obtained expressing the mode-pairs in terms 
of fundamental QNMs of a Kerr BH. We showed that it 
is the relative magnitude of the current-quadrupole mode 
S^^ and the mass-octupolc mode P^, together with the 
differences of the QNM fundamental frequencies for each 
of the dominant modes, that determine the difference 
between the recoil at the peak of the linear momentum 
flux, and the final recoil velocity, i.e., the magnitude of 
the anti-kick. 

With the final goal of improving analytic PN mod- 
els, we also explored whether simple modifications of the 
Newtonian formula for the linear momentum flux allow us 
to match the numerical results all along the binary evo- 
lution. We found that, if we treat the binary radial sep- 
aration in the Newtonian multipole modes (|37a|) - p7eP 
with an effective radius, which is computed from the nu- 
merical simulations assuming that each multipole mode 
is described by a dominant frequency (see Fig. [2]), the 
leading Newtonian modes reproduce quite well the nu- 
merical ones (see Figs. [71 [5]) up to the end of the inspiral 
phase. We also found, confirming the results in Ref. [73j . 
that a superposition of three QNMs can fit the numer- 
ical waveforms very well from the peak of the radiation 
through the RD phase. 

The tools developed in this paper will be employed to 
improve current analytic predictions for the recoil veloc- 
ity m, Elusing PN analytic models [Hi and the BOB 
approach [37l. HI, H^, HO, IH, [T^I ■ An accurate, fully ana- 
lytic description of the recoil velocity can be adopted in 
fast Monte Carlo simulations to predict recoil distribu- 
tions from BH mergers with uncertainties smaller than 
in Ref. (s^ . Those recoil distributions can in turn be 
included in simulations of hierarchical merger models of 
supermassive BHs providing more robust predictions for 
LISA. 
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FIG. 18: Flux amplitudes from tfie NEjo* run, as in Fig. g] 
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FIG. 19: Phase differences from the NEji,* run, as in Fig. 1131 



APPENDIX A: RESULTS FROM 1:4 MASS 
RATIO 
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In addition to the simulations presented in the main 
body of this paper, we have also recently analyzed a non- 
spinning system with mass ratio 1:4 {rj ~ 0.16). The 
results of this analysis are presented briefly in this ap- 
pendix, as well as in Tables Hl lIIII (labeled ap pro priately 
as NEq5^). More details can be found in Ref. |77| . 

In Fig. [18] we show the flux amplitudes from the dif- 
ferent modes, as in Fig. [5| above. We find the relative 
amplitudes almost identical to those of the NEq^^ run, 
with a slightly stronger contribution from the I'^^ mode, 
as expected from Eq. p7g[ ), which predicts a maximum 
in the amplitude for rj = 0.167. 

In Fig.[Tn|we plot the phase relations between the dif- 
ferent flux vectors, defined in Eqs. (|49ap - (|49cp . As antic- 
ipated in Sec. IVI Bl above, we find a smaller phase shift 
in ?/>'^^'^ at the transition from inspiral to ringdown for 
this more extreme mass-ratio system. The other phases 
appear to behave as expected. 
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FIG. 20: Relative contributions to the total recoil velocity 
from the different multipole mode-pairs for the NEqo* run, as 
in Fig. [13 



Lastly, in Fig. [20l we show the total recoil velocity 
along with the relative contributions from the dominant 
modes for the NEj5^ run. Again, the qualitative behav- 
ior is quite similar to the NEQij^ and NEj5^ runs, but we 
can now identify a clear trend of a smaller anti-kick for 
smaller values of t?. As mentioned above in Section fVl C) 
the amplitude of the anti-kick is most strongly dependent 
on the relative amplitudes of the 5*^^ and I^^ modes, but 
for non-spinning BH binaries, these modes both scale the 
same with mass ratio. However, the amplitude of the /^^ 
mode decreases with decreasing 77, while the amplitude 
of increases with decreasing 77, at least over the range 
considered here. Thus the amplitude of the F'^^'^^ flux in- 
creases relative to the F^^'^^ flux for more extreme mass 
ratios. From Figs. [TSl and [20l we see that the i^22,33 g^^^^ 
dominates the anti-kick, while the F'^^'^'* flux contributes 
almost nothing to it, so by increasing the relative ampli- 
tude of ^33-44^ jjg^yg effectively decreased the size of 
the anti-kick. 
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